### Alizade, Dancygier, Ditlmann 
### "National Penalties Reversed"
### Replication Code 
### Figure A6
### For questions, contact jalizade@princeton.edu

# empty environment
rm(list = ls())

#setwd("")
setwd("C:/Users/Jey/Dropbox/WZB/NaturalizationExperiment/Submission/JOP/replication_JOP/data")

# load necessary packages
library(readstata13)

# load data set
dat <- read.dta13("data_experimental.dta")

# convert treatment variable from second experiment to factor
dat$e2_treat[dat$e2_treat=="NA"] <- NA
dat$e2_treat <- factor(dat$e2_treat, levels=c("Vote Intention (Canadian)", "Control (Turkish)", "Vote Intention (Turkish)", "Integration Problems (Turkish)"))

### Figure A6a ###

# run regressions and extract coefficients
coefs_a <- lapply(paste0("pop", seq(2500, 5000, 250)), function(x) {
  mod <- lm(e1_response ~ e1_treat_turkish*dat[,x] + e1_treat_dual*dat[,x] + e1_treat_benefits*dat[,x], data=dat)
  modrob <- coeftest(mod, vcov = vcovHC(mod, "HC1"))[6:8, 1:2]
})

# create data frame with estimates
coefs_a <- do.call(rbind.data.frame, coefs_a)
names(coefs_a) <- c("coef", "se")
coefs_a$treat <- rep(c("Turkish", "Dual", "Benefits"), 11)
coefs_a$treat <- factor(coefs_a$treat, levels=c("Turkish", "Dual", "Benefits"))
coefs_a$threshold <- rep(seq(2500, 5000, 250), each=3)

# plot
plot_a <- ggplot(data=coefs_a, aes(x=threshold, y=coef))
plot_a <- plot_a + facet_wrap(~treat)
plot_a <- plot_a + geom_errorbar(aes(ymax = coef + se, ymin = coef - se), size = 1.2, width = 0)
plot_a <- plot_a + geom_errorbar(aes(ymax = coef + 1.96*se, ymin = coef - 1.96*se), size = 0.4, width = 0.2) 
plot_a <- plot_a + geom_point(aes(x = threshold, y = coef), position = position_dodge(width = 0.5), size = 3)
plot_a <- plot_a + coord_flip()
plot_a <- plot_a + scale_x_continuous(breaks=seq(2500, 5000, 250))
plot_a <- plot_a + geom_hline(yintercept = 0, linetype = "dashed")
plot_a <- plot_a + xlab("Population Threshold") + ylab(expression(paste("Interaction with Population Dummy")))
plot_a


### Figure A6b ###

# relevel treatment variable
dat$e2_treat <- relevel(dat$e2_treat, ref="Vote Intention (Turkish)")

# run regressions and extract coefficients
coefs_b <- lapply(paste0("pop", seq(2500, 5000, 250)), function(x) {
  mod <- lm(e2_response ~ e2_treat*dat[,x], data=dat)
  modrob <- coeftest(mod, vcov = vcovHC(mod, "HC1"))[6:8, 1:2]
})

# create data frame with estimates
coefs_b <- do.call(rbind.data.frame, coefs_b)
names(coefs_b) <- c("coef", "se")
coefs_b$treat <- rep(c("Vote Intention\n(Canadian)", "Control\n(Turkish)", "Integration Problems\n(Turkish)"), 11)
coefs_b$treat <- factor(coefs_b$treat, levels=c("Vote Intention\n(Canadian)", "Integration Problems\n(Turkish)", "Control\n(Turkish)"))
coefs_b$threshold <- rep(seq(2500, 5000, 250), each=3)

# plot
plot_b <- ggplot(data=coefs_b, aes(x=threshold, y=coef))
plot_b <- plot_b + facet_wrap(~treat)
plot_b <- plot_b + geom_errorbar(aes(ymax = coef + se, ymin = coef - se), size = 1.2, width = 0)
plot_b <- plot_b + geom_errorbar(aes(ymax = coef + 1.96*se, ymin = coef - 1.96*se), size = 0.4, width = 0.2) 
plot_b <- plot_b + geom_point(aes(x = threshold, y = coef), position = position_dodge(width = 0.5), size = 3)
plot_b <- plot_b + coord_flip()
plot_b <- plot_b + scale_x_continuous(breaks=seq(2500, 5000, 250))
plot_b <- plot_b + geom_hline(yintercept = 0, linetype = "dashed")
plot_b <- plot_b + xlab("Population Threshold") + ylab(expression(paste("Interaction with Population Dummy")))
plot_b
